Here, the results of analysis of the epithelial cell fraction by Seurat were visualized.
# Load libraries
library(Seurat) # version 3.0.1
library(useful) # Corner function
library(ggplot2)
library(dplyr)
library(knitr)
library(kableExtra)
library(devtools)
library(cowplot)
library(viridis)
# Load original functions
source("./Src/Visualization_for_object_of_Seurat_v3.r")myTSNEPlot(EWF, do.label = F, no.rect = F, legend_ncol = 1, label.size = 5, group.by = "stage_new", pt.size = 3,
pt.color = c("#ff4000", "#ff8c00", "#ffe500", "#b2db11", "#1b9850", "#00d3cc", "#0188a8"))# Check the batch effect by highlighting plate information
plot_grid(highlightTSNEPlot(EWF, group.by = "plate", pt.shape = 21, theme.grey = TRUE,
highlight.use = "E11_180312003"), # E11.5
highlightTSNEPlot(EWF, group.by = "plate", pt.shape = 21, theme.grey = TRUE,
highlight.use = c("E12_160809001", "E12_160809003", "E12_160809004")), # E12.0
highlightTSNEPlot(EWF, group.by = "plate", pt.shape = 21, theme.grey = TRUE,
highlight.use = c("E13_150126001", "E13_160126002", "E13_160126004")), # E13.0
highlightTSNEPlot(EWF, group.by = "plate", pt.shape = 21, theme.grey = TRUE,
highlight.use = c("E13_161012002", "E13_161012003")), # E13.5
highlightTSNEPlot(EWF, group.by = "plate", pt.shape = 21, theme.grey = TRUE,
highlight.use = c("E14_161019002", "E14_161019003")), # E14.0
highlightTSNEPlot(EWF, group.by = "plate", pt.shape = 21, theme.grey = TRUE,
highlight.use = c("E15_150126001", "E15_161014001", "E15_161014002", "E15_161014003")), # E15.0
highlightTSNEPlot(EWF, group.by = "plate", pt.shape = 21, theme.grey = TRUE,
highlight.use = c("E17_161014001", "E17_161014002", "E17_161014003")), # E17.0
ncol = 4, align = "v")df0 <- data.frame(Plate = EWF@meta.data[, "plate"], Cluster = EWF@active.ident)
df0 <- FSA::Subset(df0, df0$Plate != "kE12_180314001" | df0$Plate != "E12_160809001")
df <- df0 %>%
dplyr::group_by(Plate) %>%
dplyr::select(Cluster) %>%
table %>% as.data.frame() %>%
dplyr::group_by(Plate) %>%
dplyr::mutate(Pos = cumsum(Freq) - (Freq * 0.5))
brks <- c(0, 0.25, 0.5, 0.75, 1)
df %>%
ggplot(aes(x = Plate, y = Freq, fill = Cluster)) +
geom_bar(stat = "identity", position = "fill") +
scale_y_continuous(breaks = brks, labels = scales::percent(brks)) +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
guides(fill=guide_legend(ncol=2, byrow=TRUE))# find markers for every cluster compared to all remaining cells, report
# only the positive ones
wf.markers.res2 <- FindAllMarkers(object = EWF, only.pos = TRUE, min.pct = 0.5, thresh.use = 1, test.use = "roc", print.bar = FALSE)wf.markers.res2.top25 <- wf.markers.res2 %>% dplyr::group_by(cluster) %>% dplyr::top_n(25, myAUC) %>% dplyr::top_n(25, avg_diff) colfunc <- colorRampPalette(c("black", "#440154"))
heatmap.col <- c(colfunc(20), viridis(50))
myComplexHeatmap(EWF,
genelist = wf.markers.res2.top25,
disp.min = -2.5, disp.max = 2.5,
group.by = "labelmod.res.2_added",
heatmap.col = heatmap.col,
ann.colors = scales::hue_pal()(16),
label.set = c("Wnt9b", "Wnt7b", "Nell2", "Edar", "Wnt10b", "Shh", "Ly6d", "Krt80",
"Krt6a", "Krt1", "Aqp3", "BC100530", "Lmo1", "Tgfb2", "Nfatc1", "Fbn2",
"Adamts20", "Adamts17", "Tgfbi", "Tgfbi", "Shisa2", "Smtn", "Vdr", "Pou4f3"))genes <- list(Epithelium = c("Cdh1", "Itga6", "Krt5"),
HairGerm = c("Shh", "Lef1", "Wnt10b"),
HFSCs = c("Nfatc1", "Sox9", "Krt15"),
Infundibulum = c("Aqp3", "Klf5", "Fst"),
KeratinizedCells = c("Krt6a", "Krt1", "Krt10"),
MerkelCells = c("Krt8", "Krt20", "Sox2"))
myFP <- list()
myFP <- lapply(genes, function(i) myFeaturePlot(EWF, i, nCol = 1, title.size = 10, pch.use = 21, outline.size = 0.1,
outline.color = "grey50"))
# print(myFP[[1]])######## remove
for (i in 1:length(genes)) {
tiff(paste("Epithelium_FeaturePlot", i, ".tiff", sep = ""), width = 300, height = 1000, res=300, bg = "white")
print(myFP[[i]])
dev.off()
}gene <- list("Krt1", "Krt6a", "Ly6d", "Krt80", "Aqp3", "BC100530", "Wnt9b", "Lmo1", "Pthlh", "Nfatc1", "Tgfb2", "Fbn2",
"Adamts20", "Adamts17", "Tgfbi", "Shisa2", "Smtn", "Vdr", "Wnt10b", "Shh")
mySP <- list()
mySP <- lapply(gene, function(i) StagePlot3Mod(EWF,
stage.use = c("E11.5", "E12.0", "E13.0", "E13.5", "E14.0", "E15.0", "E17.0"),
features.use = i, title.size = 15, nCol = 7))
# print(mySP[[1]])DotPlot(object = EWF, assay = "RNA", features = c("Wnt9b", "Wnt7b", "Wnt10b", "Shh",
"Ly6d", "Krt80", "Krt6a", "Krt1", "Aqp3", "BC100530", "Tgfb2", "Nfatc1", "Fbn2",
"Adamts20", "Tgfbi", "Shisa2", "Smtn", "Vdr", "Pou4f3")) + RotatedAxis() + theme(axis.text.x = element_text(size = 16,
family = "Helvetica", face = "bold.italic"), axis.text.y = element_text(size = 16,
family = "Helvetica", face = "bold"))